In [1]:
import sympy as sp
import numpy as np
from pyNT.ode import Hamiltonian
define the potential, and plot it since we defined it...
In [2]:
A = 1.
B = 1.
C = 0.5
def Usp(x, y):
return A*sp.cos(x) + B*sp.cos(y) + C*sp.cos(x)*sp.cos(y)
def Unp(x, y):
return A*np.cos(x) + B*np.cos(y) + C*np.cos(x)*np.cos(y)
%matplotlib inline
import matplotlib
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(10, 6))
ax = fig.add_axes([.1, .1, .8, .8], projection='3d')
X, Y = np.ogrid[.0:3*np.pi:64j, .0:3*np.pi:64j]
ax.plot_wireframe(X, Y, Unp(X, Y))
Out[2]:
In [3]:
x = sp.Symbol('x')
y = sp.Symbol('y')
p = sp.Symbol('p')
q = sp.Symbol('q')
H = (p**2 + q**2)/2 + Usp(x, y)
trig_sys = Hamiltonian(
q = [x, y],
p = [p, q],
H = H)
In [4]:
traj = trig_sys.cRK(
0.01,
2**12,
np.array([[3.14, 3.14 ],
[3.14, 3.14 ],
[0.90, 0.901],
[1.61, 1.61 ]]))
In [5]:
fig = plt.figure(figsize=(9,9))
ax = fig.add_subplot(111)
ax.plot(traj[:, 0, 0], traj[:, 1, 0])
ax.plot(traj[:, 0, 1], traj[:, 1, 1])
ax.set_aspect('equal')
In [6]:
evdt0 = trig_sys.get_evdt(
h0 = 2.**(-6),
exp_range = range(4),
nsteps = 2**6,
X0 = np.array([3.14, 3.14, 0.9, 1.61]),
solver = ['Heun', 'CM2'])
evdt1 = trig_sys.get_evdt(
h0 = 2.**(-6),
exp_range = range(4),
nsteps = 2**6,
X0 = np.array([3.14, 3.14, 0.9, 1.61]),
solver = ['cRK', 'CM4'])
In [11]:
evdt2 = trig_sys.get_evdt(
h0 = 2.**(-2),
exp_range = range(4),
nsteps = 2**2,
X0 = np.array([3.14, 3.14, 0.9, 1.61]),
solver = ['CM6', 'CM8'])
In [12]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(evdt0[:, 0],
evdt0[:, 2],
marker = '.',
label = '2')
ax.plot(evdt1[:, 0],
evdt1[:, 2],
marker = '.',
label = '4')
ax.plot(evdt2[:, 0],
evdt2[:, 2],
marker = '.',
label = '6')
ax.set_xscale('log')
ax.set_yscale('log')
ax.legend(loc = 'best')
ax.grid()
Lyapunov exponent $\lambda$: \begin{equation} | \delta x(t) | \approx \exp(\lambda t) | \delta x (0) | \end{equation} One way to find the maximal Lyapunov exponent: \begin{equation} \lambda = \lim_{t \rightarrow \infty} \lim_{\delta x(0) \rightarrow 0} \frac{1}{t} \ln \frac{| \delta x(t) |}{| \delta x(0) |} \end{equation}
In [8]:
def get_data(
solver = ['Heun', 'Taylor2'],
N = 8,
h0exp = -7,
eps0 = 1e-17,
X0 = np.array([[3.14], [3.14], [0.9], [1.61]])):
data = {}
h0 = 2.**h0exp
data['x0'] = getattr(trig_sys, solver[0])(
h0, N*(2**(-h0exp)), X0)
data['x1'] = getattr(trig_sys, solver[1])(
h0, N*(2**(-h0exp)), X0)
data['eps'] = np.average(np.abs(
data['x0'][:, 0:2] -
data['x1'][:, 0:2]), axis = (1, 2))
data['t'] = h0*np.array(range(N*(2**(-h0exp))+1)).astype(np.float)
return data
In [29]:
ntraj = 128
initial_condition = np.zeros((4, ntraj), dtype = np.float)
initial_condition[0:2] = np.pi
tmp1 = np.random.randn(ntraj)
tmp2 = np.random.randn(ntraj)
initial_condition[2] = 2*tmp1 / (tmp1**2 + tmp2**2)**.5
initial_condition[3] = 2*tmp2 / (tmp1**2 + tmp2**2)**.5
In [52]:
data0 = get_data(
solver = ['Heun', 'CM2'],
N = 1024,
h0exp = -5,
X0 = initial_condition)
In [58]:
data1 = get_data(
solver = ['Heun', 'CM2'],
N = 1024,
h0exp = -6,
X0 = initial_condition)
In [59]:
fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(111)
ax.plot(data0['t'], data0['eps']/data0['t'],
label = 'CM2')
ax.plot(data1['t'], data1['eps']/data1['t'],
label = 'CM4')
ax.set_yscale('log')
ax.set_xscale('log')
ax.legend(loc = 'best')
#ax.set_ylim(0, 1)
plt.show()
In [66]:
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111)
ax.plot(data0['x0'][:, 0], data0['x0'][:, 1])
ax.set_xlim(-500, 500)
ax.set_ylim(-500, 500)
plt.show()
In [ ]: